An Analytic Solution for a Lotka-Volterra System

It is perhaps rather surprising, at least to me, that Jacobi elliptic functions can provide quite good approximate solutions to intractable systems. A previous presentation made use of the square of a Jacobi cosine function to capture the behavior of a simple yet quite nonlinear differential equation. This found application in the approximate analytic integration of orbits for power potentials with arbitrary real exponents.

This presentation continues that idea by constructing a good approximate analytic solution for the nonlinear Lotka-Volterra system

$\frac{d X}{d T} = X - X Y \frac{d Y}{d T} = - Y + X Y$

which is the restriction of a more general system to $α = 1$ . Expressions used in what follows can be found in the development of the more general system. The restricted problem has symmetries that allow a cleaner presentation of the analytic solution.

The energy integral for the restricted system is

$E = X + Y - ln X Y$

This integral allows the two dynamic variables to be expressed symmetrically in terms of each other using the Lambert W-function:

$X = - W [- \frac{e^{Y - E}}{Y}] Y = - W [- \frac{e^{X - E}}{X}]$

The maximum and minimum values for each of these variables occur when the other is equal its value at the critical point $[1, 1]$ of the system that is a center for motion in the phase plane. Given the symmetric relationships of the variables, let $Z$ represent either one as appropriate. The minima occur on the principal real branch of the Lambert W-function and the maxima on the other real branch:

$Z_{\min} = - W_{0} [- e^{1 - E}] Z_{\max} = - W_{- 1} [- e^{1 - E}]$

The period of the system can be expressed in terms of either variable. Using their common denotation, the integral expression for the period is

$T_{system} = \int_{Z_{\min}}^{Z_{\max}} \frac{d Z}{Z} [\frac{1}{[1 + W_{0} [- \frac{e^{Z - E}}{Z}]]} - \frac{1}{[1 + W_{- 1} [- \frac{e^{Z - E}}{Z}]]}]$

Since the limits of integration depend only on energy, the period is also a function of energy alone.

Knowing from this that energy is the controlling parameter of the restricted system, one can choose initial conditions of the variables to simplify the temporal dependence of the solution. One judicious was to do this is by examining combinations of the two dynamic variables. From the original pair of nonlinear differential equations, one can construct equations for the combinations of variables

$\begin{array}{l} \frac{d}{d T} (X + Y) = (X - Y) \\ \frac{d}{d T} (X - Y) = (X + Y) - 2 X Y \\ \frac{d}{d T} (X Y) = (X - Y) X Y \end{array}$

Two of the combinations, the sum and product, are easily connected through the energy integral for this restricted system, which is not possible in the more general case. It is also clear that these two combinations have extrema at the same phase-space point $X = Y$ for a given energy. That indicates one way to simplify the problem: consider initial conditions for which the two dynamic variables begin equally and calculate the energy from that single value.

As an aside, when both variables begin at unity the exact solution is identically a constant of one for both dynamic variables, while the period of the system approaches $2 π$ .

With this choice, the numerical integration of the pair of nonlinear equations over two periods along with combinations of the dynamic variables (prey is green, predators red) looks like this:

This interactive graphic is biased towards the product combination for a good reason: if the product can be approximated nicely then the differential equations become linear of first order and thus are always soluble. To wit, if one can determine some function such that $X Y = f (T)$ then the two dynamic variables are given by

$\begin{array}{l} X (T) = e^{T} [Z_{0} - \int d T e^{- T} f (T)] \\ Y (T) = e^{- T} [Z_{0} + \int d T e^{T} f (T)] \end{array}$

And indeed, a bit of playing around with Jacobi elliptic functions indicates that the product is very close to being proportional to a delta amplitude function. The motivation for this is visual inspection from having interactive graphics readily available for experimentation. If there is some deeper motivation possible, then I would love to know it myself!

It is worth noting that the sum of the two dynamic variables looks very similar to the square of a Jacobi elliptic function used in the presentation mentioned at the outset. Further playing with functions indicates that it is not difficult to find parameters for which the square of an elliptic cosine is very similar to the logarithm of a delta amplitude function. These two forms are linked via the energy integral, giving some context to the choice of fitting function.

The Jacobi delta amplitude function to be used will have the form

$dn [\frac{π}{T_{system}} \frac{K (m)}{K (0)} T, m]$

The first factor in the first argument converts the period of the system to a frequency corresponding to a circular function, remembering that the delta amplitude function has half the real period of other elliptic functions. The second factor corrects this frequency as described here to account for the presence of the elliptic parameter.

Other than the initial condition $X Y (0) = Z_{0}^{2}$ there are three other points that can be determined exactly. The two dynamic variables are equal at a second extremum of the product. The energy at this second extremum can be inverted to give

$E = 2 (Z_{ext} - ln Z_{ext}) \to Z_{ext} = - W [- e^{- E / 2}]$

which due to the symmetry in the variables corresponds to a time equal to half of the period of the system,

$X Y (\frac{T_{system}}{2}) = W [- e^{- E / 2}]^{2}$

where the principal branch of the Lambert function is used here for $Z_{0} \geq 1$ and the other real branch for $Z_{0} < 1$ , since the second extremum will be on the branch different from the common initial condition.

At the extrema of individual dynamic variables, the other variable is always unity for the restricted system. That means

$\begin{array}{l} X Y (T_{Z_{\min}}) = - W_{0} [- e^{1 - E}] \\ X Y (T_{Z_{\max}}) = - W_{- 1} [- e^{1 - E}] \end{array}$

Since the delta amplitude function is even with respect to the chosen origin, the temporal coordinates in these expressions can be taken as absolute values of the simplest possible integrals,

$\begin{array}{l} T_{Z_{\min}} = | \int_{Z_{0}}^{Z_{\min}} \frac{d Z}{Z [1 + W (- \frac{e^{Z - E}}{Z})]} | \\ T_{Z_{\max}} = | \int_{Z_{0}}^{Z_{\max}} \frac{d Z}{Z [1 + W (- \frac{e^{Z - E}}{Z})]} | \end{array}$

where now the principal branch is used for $Z_{0} \leq 1$ and the other real branch for $Z_{0} > 1$ , since the closest extrema of the individual dynamic variables remains on the same branch.

With four exact points available to match, one could concoct a variety of complicated fitting functions. In the interest of simplicity and clarity of the method, the fitting function will be taken to be

$A + B dn [\frac{π}{T_{system}} \frac{K (m)}{K (0)} T, m]$

Since the period of the system is already determined, there are only three independent parameters: the two linear coefficients and the elliptic parameter. This appears to be one of the simplest choices for a decent fit, while leaving open the possibility of more accurate form with added complexity.

The two linear coefficients can be chosen to make this function exactly match the full solution at the extrema,

$\begin{array}{c} A + B = Z_{0}^{2} \\ A + B dn [\frac{π}{2} \frac{K (m)}{K (0)}, m] = W [- e^{- E / 2}]^{2} \\ A = Z_{0}^{2} - B B = \frac{Z_{0}^{2} - W [- e^{- E / 2}]^{2}}{1 - dn [\frac{π}{2} \frac{K (m)}{K (0)}, m]} \end{array}$

with the branch chosen as described above. The two similar statements for matching at the extrema of individual variables are

$\begin{array}{l} A + B dn [\frac{π}{T_{system}} \frac{K (m)}{K (0)} T_{\min}, m] = - W_{0} [- e^{1 - E}] \\ A + B dn [\frac{π}{T_{system}} \frac{K (m)}{K (0)} T_{\max}, m] = - W_{- 1} [- e^{1 - E}] \end{array}$

can be combined with the results above for the conditions

$\begin{array}{l} \frac{Z_{0}^{2} + W_{0} [- e^{1 - E}]}{1 - dn [\frac{π}{T_{system}} \frac{K (m)}{K (0)} T_{\min}, m]} = B = \frac{Z_{0}^{2} - W [- e^{- E / 2}]^{2}}{1 - dn [\frac{π}{2} \frac{K (m)}{K (0)}, m]} \\ \frac{Z_{0}^{2} + W_{- 1} [- e^{1 - E}]}{1 - dn [\frac{π}{T_{system}} \frac{K (m)}{K (0)} T_{\max}, m]} = B = \frac{Z_{0}^{2} - W [- e^{- E / 2}]^{2}}{1 - dn [\frac{π}{2} \frac{K (m)}{K (0)}, m]} \end{array}$

Either condition can be used to produce a fit: the second, the match at the maxima of individual variables, will be chosen since the product is changing more rapidly there than at the minima. Numerical root finding will be needed in practice to determine a solution to either of these equations.

With all formulae in place, one can now visualize the fitting function for the product compared to the numerical differentiation of the pair of differential equations. For large values of the common initial condition the elliptic parameter approaches closer and closer to unity, and here one faces the possibility of inaccuracy using a machine-precision numerical library. The same is true for small values of the common initial condition, where the elliptic parameter becomes a very large negative number. To avoid such inaccuracies the common initial condition will be restricted to a range where the numerical solutions are accurate.

The fitting function compared to the numerical differentiation of the product looks like this:

The fit is not exact, but one can see by switching to the difference between it and the accurate numerical solution that it is quite good. Really not bad at all for such a relatively simple fitting function.

Another way to test the basic fit is to construct a differential equation for the product combination

$\begin{array}{c} \frac{d^{2}}{d T^{2}} ln (X Y) = \frac{d}{d T} (X - Y) = (X + Y) - 2 X Y \\ \frac{d^{2}}{d T^{2}} ln (X Y) = E + ln (X Y) - 2 X Y \end{array}$

and visualize how closely the two sides match. With the abbreviation $ω = \frac{π}{T_{system}} \frac{K (m)}{K (0)}$ , the second derivative of the logarithm of the fitting function is

$\begin{array}{l} \frac{d^{2}}{d T^{2}} ln [A + B dn (ω T, m)] = \frac{d}{d T} \frac{ω B {dn}^{'} (ω T, m)}{A + B dn (ω T, m)} \\ = \frac{ω^{2} B {dn}^{''} (ω T, m)}{A + B dn (ω T, m)} - \frac{ω^{2} B^{2} [{dn}^{'} (ω T, m)]^{2}}{[A + B dn (ω T, m)]^{2}} \\ = \frac{ω^{2} B [(2 - m) dn (ω T, m) - 2 {dn}^{3} (ω T, m)]}{A + B dn (ω T, m)} \\ - \frac{ω^{2} B^{2} m^{2} {sn}^{2} (ω T, m) {cn}^{2} (ω T, m)}{[A + B dn (ω T, m)]^{2}} \end{array}$

This expression could be simplified further by writing it entirely in terms of the delta amplitude function, but for a numerical evaluation this is unnecessary. With right- and left-hand sides colored as starboard and port, the two look like this:

The two functions are of course not equal because the analytic solution is merely approximate. That they are generally about the same is an indication that the approximation has merit.

One thing that remains is to apply the fitting function for the product of dynamic variables to evaluate each individual one using the integral forms already given. For the predators this is straightforward and produces an expected result over a single period:

To prevent the integration from slowing down the browser excessively, the approximation is evaluated at a reduced number of selected points of the numerical integration and cached in incremental steps.

The application of the approximation to the prey population over a single period is more problematic:

Virtually the same code is used for this interactive graphic, with adjustments required to evaluate the second dynamic variable, yet the outcome is clearly inaccurate. Since all quantities involved are well above tolerance values used in the integration, the result is quite puzzling. It is likely due to an inherent numerical stability in how this particular variable is calculated, but remains an open problem to be resolved.

There are a couple ways around this predicament. Since the predator population variable is accurate, the prey population variable can be evaluated by the simple relation

$X (T) = \frac{X Y (T)}{Y (T)}$

replacing the product with the fitting function already determined. One can also use the exact symmetry

$X (T) = Y (- T) = Y (T_{system} - T)$

for this system to evaluate the prey population variable.

It is also worth nothing that one can rearrange differential equations above for the exact results

$\begin{array}{l} X (T) = \frac{1}{2} [\frac{d^{2}}{d T^{2}} ln (X Y) + \frac{d}{d T} ln (X Y) + 2 X Y] \\ Y (T) = \frac{1}{2} [\frac{d^{2}}{d T^{2}} ln (X Y) - \frac{d}{d T} ln (X Y) + 2 X Y] \end{array}$

Plotting these combinations using derivatives above of the fitting function produces graphics that follow the general form of each individual dynamic variable qualitatively but not very accurately. Apparently both integration and differentiation of the approximate fitting function can amplify discrepancies between the numerical integration of the differential equations and the analytic approximation.